official_birth_records = read_feather(path = str_c(IO$out_Rdata,"official_birth_records.feather"))official_birth_records = official_birth_records %>% mutate(births_thousands = births/1000)
official_births_summary = official_birth_records %>%
group_by(country_area, country_area_col) %>%
summarize(min_births = min(births_thousands),max_births = max(births_thousands),mean_births = mean(births_thousands),
max_date = max(date))
ref_date = as.Date("2020-01-01")
g_births = ggplot(official_birth_records, aes(x = date, y = births_thousands, col = country_area_col))
g_births = g_births +
#### reference for variations
geom_segment(data = official_births_summary, aes(x = ref_date, xend = ref_date, y = mean_births, yend = mean_births-5*mean_births/100),
size = 1.5, alpha = 0.7)+
# geom_text(data = official_births_summary, aes(x = ref_date, y = min_births-5*mean_births/100), label = "5% change from the mean",
# hjust = 0, vjust = 0, nudge_x = 100)+
#### uncorrected births
geom_line(aes(y = births_original_numbers/1000), col = "gray80")+
#### corrected births
geom_line()+ #geom_point(size = 0.5)+
#### COUNTRIES
geom_text(data = official_births_summary, aes(x = ref_date, y = max_births *1.05, label = country_area, fontface = 2), vjust = 0, hjust = 1)+
#### settings
scale_color_identity()+
scale_x_date(date_minor_breaks = "1 year", date_labels = "%Y")+
scale_y_continuous(expand = expansion(mult = c(0,0.2)))+ # sec.axis = sec_axis(~ (. - max(.))/max(.)*100, breaks = seq(-30,30,by = 10), name = "% change from the max")
ylab("Births (monthly in thousands)")+ xlab("Date")+
facet_grid(country_area ~., scale = "free")+
theme(strip.text.y = element_blank())
g_births### change
users = read_feather(path = paste0(IO$output_clue, "users.feather"))
n_users_per_country = users %>% group_by(country_area) %>% summarise(n_users = n()) %>%
mutate(country_area_col = dict$country_area$country_area_col[match(country_area, dict$country_area$country_area)],
country_area = factor(country_area, levels = dict$country_area$country_area))
g_n_users = ggplot(n_users_per_country, aes(x = country_area, y = n_users/1000, fill = country_area_col))+
coord_flip()+
geom_bar(stat = "identity")+
scale_fill_identity()+
xlab("")+ylab("K# of app users")+
theme(strip.text = element_blank(),
axis.text.y = element_blank())+
facet_grid(country_area ~ ., scale = "free")
g_n_userspp = read_feather(path = str_c(IO$output_clue,"pop_indicators/relative_sexual_indicators_detrended.feather"))
pp_long = pp %>% filter(BC == "all") %>%
dplyr::select(date, country_area, n_users, n_sex, n_prot_sex, n_unprot_sex) %>%
rename(total_sex = n_sex, protected_sex = n_prot_sex, unprotected_sex = n_unprot_sex) %>%
pivot_longer(cols = c("total_sex","protected_sex","unprotected_sex"), names_to = "sex_type", values_to = "n_log") %>%
mutate(sex_type_str = sex_type %>% str_replace(.,"_"," ") %>% str_to_title())
pp_sum = pp_long %>%
group_by(date, sex_type, sex_type_str) %>%
summarize(n_users = sum(n_users),
n_log = sum(n_log)) %>%
mutate(rel_counts = n_log / n_users) %>%
arrange(sex_type , date) %>% ungroup() %>% group_by(sex_type) %>%
mutate(
t = row_number(),
rel_trend = loess(rel_counts ~ t) %>% predict(),
sex_type_str = sex_type_str %>% factor(., levels = c("Total Sex","Protected Sex","Unprotected Sex")))
g_sex_freq = ggplot(pp_sum, aes(x = date, y = rel_counts, col = sex_type_str))+
geom_line()+
geom_line(aes(y = rel_trend))+
ylab("Sex frequency (daily)")+xlab("Date")+
expand_limits(y = 0)+
scale_color_manual(name = "",values = c("gray40","aquamarine3","indianred1"))+
theme(legend.position = "top",
legend.direction = "horizontal",
legend.title = element_blank(),
legend.background = element_rect(fill = "white", color = "transparent", size = 0.5)) #color = "gray90",
g_sex_freqpp_sum = pp_sum %>%
group_by(sex_type) %>%
mutate(
mean_rel_counts = mean(rel_counts),
detrended_counts = rel_counts/rel_trend * mean_rel_counts)
g_rel_sex = ggplot(pp_sum, aes(x = date, y = detrended_counts, col = sex_type_str))+
geom_line()+
ylab("Detrended sex frequency (daily)")+xlab("")+
expand_limits(y = 0)+
scale_color_manual(name = "",values = c("gray40","aquamarine3","indianred1"))+
theme(legend.position = "top",
legend.direction = "horizontal",
legend.title = element_blank(),
legend.background = element_rect(fill = "white",color = "transparent", size = 0.5)) # gray90
g_rel_sex# import png
clue_screen = image_read("../Figures Tables Media/Media/Clue_app_screens.png")
clue_screen_ratio = image_info(clue_screen)$height / image_info(clue_screen)$width
padding = 20
g_clue_screen = ggplot()+ coord_fixed(ratio = clue_screen_ratio)+
background_image(clue_screen)+theme(plot.margin = margin(t=padding, l=padding, r=padding, b=padding, unit = "pt"))fig_1 = arrangeGrob(
g_births,
g_clue_screen,
g_rel_sex,
ncol = 2, nrow = 2,
widths = c(1.5, 1),
heights = c(0.55,0.45),
layout_matrix = cbind(c(1,1), c(3,4)))
fig_1 = as_ggplot(fig_1) +
draw_plot_label(label = letters[1:3], size = 15,
x = c(0, 0.6, 0.6), y = c(1, 1, 0.45))
fig_1# load Figure 2 data
load(file = str_c(IO$p_outputs,"sex_models.Rdata"), verbose = TRUE)## Loading objects:
## sex_models
## sex_models_df
plotlist = list()
n_items = c()
cas = unique(sex_models_df$country_area)
ok = foreach(ca = cas) %do% {
cat(ca, "\n")
# retrieving the model for this country and only plotting for BC all and sex_type all
j = which((sex_models_df$country_area == ca) & (sex_models_df$BC == "all") & (sex_models_df$sex_type == "unprot_sex"))
glm_sex_behavior = sex_models[[j]]$model
# Visualisation of the coefficient
g_coef = ggplot_sex_activity_glm_coefficient(model = glm_sex_behavior, show_intercept = FALSE, show_weekly_patterns = FALSE, ylim = c(-0.4, 1))
g_coef = g_coef +
ggtitle(ca)+
theme(axis.text.x = element_text(size = 0.8))
print(g_coef)
plotlist[[ca]] = g_coef
n_items = c(n_items, nrow(g_coef$data))
}## Brazil - Central-West
## Brazil - Northeast
## France
## United Kingdom
## United States - California
## United States - Northeast
fig_2 = ggarrange(
ggarrange(
plotlist[[which(cas == "Brazil - Central-West")]],
plotlist[[which(cas == "United States - Northeast")]],
labels = c("a","b"),
ncol = 2, nrow = 1,
widths = n_items[which(cas %in% c("Brazil - Central-West","United States - Northeast"))] /
sum(n_items[which(cas %in% c("Brazil - Central-West","United States - Northeast"))])
),
ggarrange(
plotlist[[which(cas == "France")]],
plotlist[[which(cas == "United Kingdom")]],
labels = c("c","d"),
ncol = 2, nrow = 1,
widths = n_items[which(cas %in% c("France","United Kingdom"))] /
sum(n_items[which(cas %in% c("France","United Kingdom"))])
),
nrow = 2, ncol = 1
)
fig_2G_table = read_feather(path = str_c(IO$out_Rdata,"Gestation_par_table.feather"))# t = seq(0,4*pi, by = pi/20)
# f = sin(t)
# plot(t, f, type = "l")
t = seq(30*7, 45*7)
d = foreach(ca = G_table$country_area, .combine = bind_rows) %do% {
this_ca = G_table %>% filter(country_area == ca)
res = this_ca[rep(1,length(t)),]
res = res %>% mutate(t = t, d = dnorm(t, mean = this_ca$G*7, sd = Gsd0), G = this_ca$G*7)
res
}
br = seq(min(t), max(t), by = 14)
d = d %>%
mutate(country_area_wrapped = str_replace(country_area, " - ","\n"),
country_area_col = dict$country_area$country_area_col[match(country_area, dict$country_area$country_area)])
g_d = ggplot(d, aes(x = t, y = d, fill = country_area_col))+
geom_area(alpha = 0.6)+
geom_vline(aes(xintercept = G, col = country_area_col))+
scale_color_identity()+
scale_fill_identity()+
scale_x_continuous(breaks = br , labels = br/7)+
scale_y_continuous(breaks = NULL)+
xlab("G, gestation duration\n(weeks)")+
ylab(expression(paste("Density, d",tau)))+ #expression(paste("Value is ", sigma,",", R^{2},'=0.6'))
facet_grid(country_area_wrapped ~ ., scale = "free")+
theme(strip.text.y = element_text(angle = 0, hjust = 0))
g_dggplot(d, aes(x = t, y = country_area, alpha = d))+
geom_tile()+
scale_y_discrete(limits = rev(levels(G_table$country_area)))+
xlab("G, gestation duration\n(weeks)")+
ylab("Density")+
guides(alpha = FALSE)+
scale_alpha_continuous(range = c(0,1))+
scale_x_continuous(breaks = br , labels = br/7)model = image_read("../Figures Tables Media/Media/Figure 3-01.png")
model_ratio = image_info(model)$height / image_info(model)$width
g_model = ggplot()+ coord_fixed(ratio = model_ratio)+
background_image(model)
models = image_read("../Figures Tables Media/Media/Figure 3-02.png")
models_ratio = image_info(models)$height / image_info(models)$width
g_models = ggplot()+ coord_fixed(ratio = models_ratio)+
background_image(models)fig_3 = ggarrange(
ggarrange(
g_model,
g_d,
labels = c("a","b"),
ncol = 2, nrow = 1,
widths = c(1,1)
),
ggarrange(
g_models,
labels = c("c"),
ncol = 1, nrow = 1,
widths = 1
),
nrow = 2, ncol = 1,
heights = c(1,1)
)
fig_3births = read_feather(path = str_c(IO$out_Rdata, "simulated_births.feather"))
births_STL = read_feather(path = str_c(IO$out_Rdata, "simulated_births_seasonal_trends.feather"))
opt_par_df = read_feather(path = str_c(IO$out_Rdata, "optimal_parameters_and_AIC.feather"))
ca_levels = dict$country_area$country_area
births = births %>% mutate(country_area = factor(country_area, levels = ca_levels))
births_STL = births_STL %>% mutate(country_area = factor(country_area, levels = ca_levels))
opt_par_df = opt_par_df %>% mutate(country_area = factor(country_area, levels = ca_levels))births = births %>%
mutate(country_area_wrapped = str_replace(country_area," - ","\n") %>%
factor(., levels = dict$country_area$country_area %>% str_replace(.," - ","\n")))
g_ts = ggplot(births %>% filter(sex_type == "unprot_sex", BC == "all"), aes(x = year_month, y = sim_births, col = model))
g_ts = g_ts +
geom_line(aes(y = births), col = "black", size = 1)+
geom_line(size = 0.5)+
guides(col = FALSE)+
xlab("Date")+ylab("")+
facet_grid(country_area_wrapped ~ model, scale = "free" ,switch = "y")+
theme(strip.placement = "outside",
strip.text.y.left = element_text(angle = 0, hjust = 0.5),
strip.background.y = element_rect(fill = "gray90", color = NA))+
ggtitle("Measured vs. simulated births")
g_tsg_st = ggplot(births_STL, aes(x = month, y = sim_births_seasonal, col = model))
g_st = g_st +
geom_hline(yintercept = 0, col = "gray80")+
geom_line(aes(y = births_seasonal), col = "black", size = 1)+
geom_line()+
scale_x_continuous(breaks = seq(0,12,by = 3))+
guides(col = FALSE)+
xlab("Calendar months")+ylab("Seasonal trends")+
facet_grid(country_area ~ model, scale = "free_y")+
theme(strip.text.y = element_blank())+
ggtitle("Seasonal trends")
g_stopt_par_df = opt_par_df %>% group_by(country_area, BC, sex_type) %>%
mutate(AIC_best_model = min(AIC),
dAIC = AIC - AIC_best_model)
g_aic = ggplot(opt_par_df %>% filter(BC == "all", sex_type == "unprot_sex"),
aes(x = model, y = country_area, fill = dAIC))
g_aic = g_aic +
geom_tile()+
scale_fill_gradient(name = "", low = hsv(0.58,0.4,1), high = hsv(0.02,0.9,0.9))+
xlab("")+ylab("")+
#guides(fill = FALSE)+
facet_grid(country_area ~ model, scale = "free")+
ggtitle("Diff. with best AIC")+
theme(axis.text = element_blank(),
strip.text.y = element_blank())
g_aicg_best_aic = ggplot(opt_par_df %>% filter(BC == "all", sex_type == "unprot_sex", model == "A"),
aes(x = country_area, y = AIC_best_model))
g_best_aic = g_best_aic +
geom_bar(stat = "identity")+ coord_flip()+
xlab("")+
facet_grid(country_area ~ model, scale = "free")+
theme(axis.text = element_blank(), strip.text.x = element_blank())
g_best_aicopt_par_df = opt_par_df %>%
mutate(country_area_col = country_area_col %>% factor(., levels = dict$country_area$country_area_col))
g_F = ggplot(opt_par_df %>% filter(model != "A", sex_type == "unprot_sex", BC == "all"),
aes(x = Tp, y = alpha, col = country_area_col))
g_F = g_F +
geom_segment(aes(xend = Tp, yend = 0, linetype = model))+
geom_point()+
scale_color_identity(guide = "legend",
labels = dict$country_area$country_area,
name = "countries/areas" )+
scale_x_continuous(limits = c(0,1),breaks = seq(0,3/4,by = 1/4), labels = c("Jan","Apr","Jul","Oct"))+
ylab("Relative amplitude")+ xlab("Fertility peak time")+
coord_polar(theta = "x", start = 0)
g_Fopt_par_df = opt_par_df %>%
mutate(
country = str_split_fixed(country_area, " - ",2)[,1],
Tp_seas = (Tp-10/365 - 0.5*(country == "Brazil"))%%1)
g_F_seas = ggplot(opt_par_df %>% filter(model != "A", sex_type == "unprot_sex", BC == "all"),
aes(x = Tp_seas, y = alpha, col = country_area_col))
g_F_seas = g_F_seas +
geom_segment(aes(xend = Tp_seas, yend = 0, linetype = model))+
geom_point()+
scale_color_identity(guide = "legend",
labels = dict$country_area$country_area,
name = "countries/areas" )+
scale_x_continuous(limits = c(0,1),breaks = seq(0,3/4,by = 1/4), labels = c("Winter\nSolstice","Spring\nEquinox","Summer\nSolstice","Autumn\nEquinox"))+
ylab("Relative amplitude")+ xlab("Fertility peak time")+
coord_polar(theta = "x", start = 0)
g_F_seasfig_4 = ggarrange(
ggarrange(
g_ts,
g_st,
g_aic,
labels = letters[1:3],
ncol = 3, nrow = 1,
widths = c(3,1,0.7)
),
ggarrange(
g_F,
g_F_seas,
ggplot(),
labels = c("e","f",""),
ncol = 3, nrow = 1,
widths = c(1,1,1), heights = 1
),
nrow = 2, ncol = 1,
heights = c(2,1)
)
fig_4Figures_path = "../Figures Tables Media/Figures/"
scale = 2
ggsave(plot = fig_1, filename = str_c(Figures_path, "F1.pdf"),
width = 17.5, height = 12, units = "cm", scale = scale)
ggsave(plot = fig_2, filename = str_c(Figures_path, "F2.pdf"),
width = 17.5, height = 10, units = "cm", scale = scale)
ggsave(plot = fig_3, filename = str_c(Figures_path, "F3.pdf"),
width = 8, height = 7, units = "cm", scale = scale)
ggsave(plot = fig_4, filename = str_c(Figures_path, "F4.pdf"),
width = 17.5, height = 13, units = "cm", scale = scale)